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SYMBOLS 


I 

Some  of  the  notation  is  borrowed  from  References  1 and  2 , 
where  more  detailed  descriptions  may  be  found. 

m 

[A]  ~ material  stiffness  defined  by  Eq.  (31) 

~ tensor  defined  by  Eq,  (79) 

[B]  ~ defined  by  Eq.  (8) 

B^j  - left  Cauchy-Green  tensor 

D-O-P  ~ degree-of-polymerization 

e^  ~ defined  by  Eq.  (16) 

E ~ Green -St,  Venant  strain  tensor 

f -■  body  force 

m 

{f}  ~ element  force  vector  defined  by  Eq.  (19) 

F - deformation  gradient 

{F}  - global  force  vector  defined  by  Eq.  (24) 

g ~ metric  tensor 

G ~ scaler  shear  modulus 

m . • 

[G]  - incidence  array  for  the  m^n  element,  defined  by  Eq.  (15) 

IB  tr  B^j 

V 


1 


J 


Det  F 


[k] 

K 

[K] 

K“Y6fi 


[N] 

{q> 

{Q} 


{T} 

T 


u 

V.W. 


value  of  J corresponding  to  zero  stress 
element  stiffness  defined  by  Eq.  (33) 
scaler  bulk  modulus 

global  stiffness  matrix  defined  by  Eq,  (34) 
elasticity  tensor 
defined  by  Eq.  (7) 

current  coordinates  of  the  nodes  of  the  m^^  element 
global  vector  representing  current  node  coordinates 
surface  traction 
time 

element  force  vector  defined  by  Eq.  (19) 
global  force  vector  defined  by  Eq.  (24) 

Cauchy's  stress  tensor 

first  Piola-Kirchhoff  stress  tensor 

second  Piola-Kirchhoff  stress  tensor 

displacement 

virtual  work 


Vi 


I.  INTRODUCTION  AND  SUMMARY 


Two  studies  are  presented  in  this  report.  The  first  (pre- 
sented in  Section  II)  is  the  development  of  a new  finite  element 
solution,  and  the  second  is  an  analysis  of  residual  stress  in 
filamentary  composite  material. 

The  new  feature  of  the  finite  element  solution  is  in  the 
formulation  phase  of  the  development;  existing  procedures  are 
used  for  the  numerical  solution.  The  technique  presented  is 
restricted  to  problems  usually  described  as  quasi-static  or 
creeping  flow,  but  it  is  not  restricted  by  small  strain  or  rota- 
tion, nor  is  it  specialized  to  a particular  theory  of  consti- 
tutive behavior.  The  finite  element  equilibrium  equations  are 
derived  from  the  principle  of  virtual  work  — specifically  the 
representation  which  uses  the  deformation  gradient  as  the  measure- 
of-strain  which  makes  this  technique  different  from  existing 
techniques . One  advantage  of  the  deformation  gradient  is  that  it 
is  expressible  as  a linear  function  of  the  node  coordinates 
exactly  for  a simplex  element.  Another  is  that  the  deformation 
gradient  is  an  adept  measure-of-strain  for  both  solid  and  fluid 
behavior.  Therefore,  the  new  technique  may  offer  advantages  for 
the  study  of  solid-fluid  interactions  or  for  material  behavior 
which  fits  neither  of  the  concepts  of  solids  or  fluids. 

Following  the  development  of  the  technique,  four  sample 
problems  are  studied.  A comparison  of  numerical  predictions  with 
some  exact  analytic  solutions  is  presented,  and  the  comparison 
is  excellent. 

Quite  a different  subject  is  addressed  in  Section  III: 
residual  stress.  Its  development  in  polymeric  filamentary  rein- 
forced composite  materials  is  discussed  in  terms  of  differential 
dimensional  changes  that  occur  during  processing.  Two  phenomena 
which  contribute  to  the  development  of  residual  stress  are  con- 
sidered: dimensional  change  of  the  resin  due  to  polymerization 
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and  differential  thermal  expansion  of  the  resin  and  reinforcement 
upon  cooling  from  the  cure  temperature  to  room  temperature.  In 
the  theory  of  composites  in  use  presently  the  effects  of  polymer- 
ization are  neglected;  it's  a goal  of  the  study  presented  here  to 
assess  this  assumption.  This  assumption  is  not  supported  by  this 
study — neither  on  the  basis  of  available  experimental  studies  of 
residual  stress  nor  on  the  basis  of  numerical  calculations  which 
are  presented  in  this  report. 
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II . FINITE  ELEMENT 

1 . Development  of  Formulation  and  Numerical  Solution  Technique 

It  is  convenient  to  start  with  the  definition  of  Virtual 
Work  for  a continuvim^  : 

V.W.  » (p  S^6z^  da  + 

•'s 

The  Principle  of  Virtual  Work  (for  quasistatic  motion  or  creeping 
flow)  requires 

V.W  = 0 

The  form  represented  by  Eq.  (1)  encompasses  arbitrarily  large 
displacements  and  all  models  of  constitutive  behavior.  Unfor- 
tunately, it  is  referred  to  the  current  configuration  (Eulerian 
description)  completely,  which  is  not  convenient  for  solids,  but 
which  is  used  for  fluids.  Generally,  one  uses  a form  of  V.W. 
which  is  referred  to  a reference  configuration  for  solids . The 
appropriate  quantities  are  defined  in  terms  of  their  transfor- 
mations from  the  above  quantities  and  then  Eq.  (1)  transforms 
to  the  following; 


dv 


(1) 


V.W.  = ^ da^r  + J dV  (2) 

where  6Ejjj^  = ^^k;m) 

This  is  the  so-called  Lagrangian  description,  and  it  represents 
a logical  starting  place  for  the  development  of  generally  all  of 
the  published  developments  for  solids,  at  least  for  the  so-called 
direct  stiffness  models. 
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Another  form  for  V.W,  in  terms  of  a reference  configuration 
can  be  developed.  In  this  case  the  first  Piola-Kirchhoff  repre- 
sentation of  the  stress  tensor  (Tj^)  is  used  rather  than  the 
second  (T)  . As  for  the  previous  form,  the  third  form  is 
developed  by  substituting  the  definitions  of  the  appropriate 
quantities  into  Eq.  (1): 


V.W.  = ^ 


T,kK, 


dA^  + 


/.[■ 


dV 


(3) 


It  is  this  form  of  virtual  work  which  is  used  in  the  present 
development . 

Now  the  kinematics  of  an  element  will  be  described.  Consider 
a single  tetrahedral  element  subject  to  homogeneous  deformation 
(i.e.,  measured  relative  to  the  reference  configuration).  This 
deformation  is  defined  as  follows : 

2^  - ei  + fJz’'  (4) 

where  the  e^  's  and  ' s are  constant  over  the  element  vol- 
\ame.  The  coordinates  q^  are  identified  with  the  current  coordi- 
nates of  each  node  where  the  superindex  "n"  refers  to  one 

of  the  nodes : 


Therefore , 


q^  = for  n = 1,  2,  3,  A 

evaluating  Eq.  (4)  at  the  nodes  gives 
qi  = for  n = 1,  2,  3,4 


(5) 


(6) 


Solving  Eq.  (6)  for  e^  ' s and  Fi  ' s , it  follows  that  they 

n • 

are  linear  functions  of  q^  ' s whose  coefficients  depend  on  the 
nodal  positions  in  the  reference  configuration  only.  After  these 
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solutions  are  substituted  into  Eq.  (4),  it  follows  that  the  cur- 
rent coordinate  of  a material  point  within  the  element  depends 
linearly  on  both  its  material  coordinate  and  the  current  positions 
of  the  nodal  points  ' s^  . Eq . (4)  may  be  presented  in  the 

following  convenient  form: 


{z}  = Sz  f = [N]  {q} 

z3 

where  N is  a 3 x 12  matrix  and  q is  a 12-element  vector. 
Also , the  elements  of  N are  linear  functions  of  the  material 
coordinates  (Z^  ' s)  and  are  independent  of  q^  ' s , Now  the 
deformation  gradient  is  computed  from  Eq.  (7)  and  can  be  repre- 
sented in  the  following  matrix  form: 


{F}  = 


= i 


[B]  {q} 


where  B is  a 9 x 12  matrix  whose  elements  are  independent  of 
both  the  material  coordinates  ' s and  qi  ' s . 

Now  consider  an  arbitrary  variation  of  the  current  positions 
of  the  nodal  points  ^q^  ' s^  and  derive  the  associated  variations 
of  the  current  coordinates  and  the  deformation  gradient  from 
Eqs.  (7)  and  (8): 

{6z}  = [N]  {6q}  (9) 

{6F}=[B]{6q}  (10) 

This  simple  result  follows  the  fact  that  N and  B are  inde- 
pendent of  the  current  configuration. 


It  is  convenient  to  represent  Eq.  (3)  in  the  following  matrix 


form: 


v.w. 


(11) 


£{6z)T[Tf  dA^]  + 


{6z}^Po{f}  - {«F) 


dV 


Now  introduce  the  homogeneous  deformation  which  is  represented  by 
Eqs.  (9)  and  (10): 


V.W. 


{6q} 


{[f  dA,}  . J 


[N]^p, 


{f}  - [B] 


'(tr)  'J’] 


(12) 


This  formula  for  the  V.W.  of  an  element  subjected  to  homogeneous 
deformation  is  otherwise  completely  general  and  exact;  no  approx- 
imations have  been  introduced  which  would  lim.it  arbitrarily  large 
displacements  and  strains  and  the  constitutive  behavior.  V.W. 
must  be  zero  for  arbitrary  variations,  therefore  the  bracketed 
term  must  be  zero. 


CN]^Pj^{f} 


dV  + 


(13) 


Once  the  constitutive  equation  is  defined,  which  relates  the 
stress  tensor  (Tr)  to  the  current  coordinates  of  the  nodes  (q)  , 
the  above  equation  yields  the  equilibrium  equation  for  an  element. 
Fortunately,  there  have  been  many  recent  contributions  to  the 
field  of  nonlinear  constitutive  equations , including  the  formu- 
lation of  general  invariance  requirements , which  are  useful  in  the 
development  of  the  formulation  required  here. 


Very  generally,  the  constitutive  equation  for  a simple 
material  may  be  represented  in  terms  of  a functional  (H)  ^ : 


,kM 

R 


HkM 


pi(t) 


(14) 


It  is  useful  to  discuss  the  constitutive  equation  in  these  general 
terms  so  that  the  discussion  applies  to  any  ideal  model.  Of  course. 
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the  specifics  of  the  equations  sketched  here  can  be  developed  only 
after  an  ideal  model  is  selected.  Generally,  an  anproximate 
representation  of  Eq.  (14)  will  be  required  to  develop  a solution 
to  Eq.  (13)  iteratively.  In  the  case  of  elastic  behavior  "exact" 
solutions  can  be  developed  by  the  method  of  successive  approxi- 
mations even  for  moderately  large  increments.  This  is  because 
the  stress  is  a function  of  the  deformation  gradient  and  does  not 
depend  on  the  path.  For  other  ideal  models,  e.g.,  plastic  and 
viscoelastic,  the  stress  does  depend  on  the  path.  For  these 
cases  n\americal  solutions  are  developed  using  small  incremental 
steps  in  time  (i.e.,  marching  techniques).  With  this  backgroimd 
the  solution  technique  used  in  this  development  is  a combination; 
incremental  time  steps  are  taken  and  the  solution  is  developed  at 
the  end  of  each  step  by  successive  approximations.  Iterative 
techniques  that  have  been  used  for  finite  element  solutions  are 
described  by  Stricklin  et  al.*.  The  modified  Newton-Raphson 
method  has  been  used  with  considerable  success  and  is  used  for  the 
present  development. 

Now  consider  a finite  body  (B)  of  volume  V and  surface 
S . Partition  B into  m tetrahedral  elements  and  denote  the 
current  coordinates  of  the  nodes  by  {Q}  . If  B has  N degrees - 
of-freedom,  it  follows  that  {Q}  is  an  N-vector  . 

The  connectivity  of  the  model  is  established  by  the  mapping 

{q}  » [G]  {Q}  (15) 

where  {q}  represents  the  current  nodal  coordinates  for  the  m^h 

element,  [S]  is  the  incidence  array  for  the  m^^  element  defined 
m 

by  Oden** , and  [G]  is  a six  by  N matrix  for  the  present  case. 
Also,  since  the  element  coordinates  and  the  global  coordinates  are 
referred  to  the  same  Cartesian  coordinate  system,  it  follows  that 


m 

each  row  of  [G]  has  one  unit  element  and  the  others  are  zero. 
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Now  consider  the  global  equllibrltm  equations.  It's  assumed 
that  the  deformation  Is  piece-wise  linear  on  B . The  deformation 
of  the  m^^  element  Is  represented  as  follows : 


If 

+ 


(16) 


It  follows  Eqs.  (7)  and  (8)  that  for  the  m^^  element 


m m 
{2}  - [N]  {q} 

(17) 

ID  ID  ID 

{F}  = [B]  {q} 


-ID  m ID 

where  {q}  represents  the  current  coordinates,  and  [N],CB] 
represent  CN],[B]  for  the  m^^  element. 


The  V.W.  represented  by  Eq.  (3)  must  be  zero.  After  sub- 
stituting Eq.  (17)  Into  Eq.  (3)  the  V.W.  of  B may  be  expressed 
as  follows : 


V.W. 


(18) 


Tf 

where  Sm  Is  the  surface  of  the  m^“  element  which  Is  part  cf  S 
of  B . It's  convenient  to  define  element  forces: 


(19) 
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Then  Eq.  (18)  has  the  following  representation: 


V.W.  - 2 {6?}^  - {t}JJ  (20) 

Now  consider  a virtual  displacement  of  the  global  coordinates 
{Q}  . It  follows  Eq.  (15)  that 

{6q}  = [G]  {6Q}  (21) 

Eq.  (21)  is  substituted  into  Eq.  (20). 

M 

T n T 

V.W.  - 

Now  require  that  V.W.  be  zero  for  arbitrary  variations  {6Q}  , it 
follows  Eq.  (22)  that 

{T}  » {F}  (23) 

where 

M 

m T ® 

{T}  5 2^[G]^{t} 

(24) 

M 

m T m 

{F}  E 2,CG]^{f} 
in*l 


Eq.  (23)  represents  the  global  equilibrium  equations. 


9 


For  creeping  flow  all  of  the  preceding  equations , which  are 
functions  of  time,  must  be  satisfied  for  any  time.  If  the  indi- 
vidual quantities  are  evaluated  at  two  different  times , say  t^ 
and  . then  increments  may  be  defined  as  follows : 


m r™/  x'l 
{Aq}  = 

(iQ)  E 


(25) 


n in  in 

Since  [N] , [B] , and  [G]  are  independent  of  the  current  coordi- 


nates, they  are  constant  (i.e.,  independent  of  time).  It  follows 
Eq.  (25)  and  previous  equations  that 


{Lz}  - [N]  {Aq} 


m 

{AF} 


m m 
CB]  {Aq} 


m 

{At}  = 

f 

t> 

»-h9 

II 

^m 

d V. 


m 


LR 


■'V. 


[ni’^PrCa?}  dVjj, 


(26) 


m 


M 


{AT} 


in  T ® 
2^  [G]^{At} 
m-l 


M 


{AF} 


m fp  in 

2^  [Gr{Af} 

in»l 


and  the  incremental  equilibrium  equations  follow  Eq.  (23). 


{AT}  - (AF) 


(27) 


Eq.  (27)  represents  a system  of  N nonlinear  equations.  The 
vector  {AF}  , which  represents  an  increment  of  applied  loading, 
is  specified.  It  follows  Eqs . (26)**®,  (17)*,  (15),  and  (14)  that 
the  N-vector  {AT}  depends  on  the  history  of  {Q}  . 

To  solve  the  incremental  solution  using  the  modified  Newton- 
Raphson  method  the  following  recurrence  formula  is  used. 


r+1 
{ AQ  } 


where  the  estimate  of  {AQ}  after  the  r^^  iteration  is 


and 


R 

R V'  r 

{AQ}  E 

r-1 


R 

{AT} 


|AT(Q(ti)  + AQ)j 


(28) 


(29) 


[K(ti)]  is  an  invertible  matrix  which  represents  the  stiffness 
of  the  global  finite  element  model  at  t = t^  . For  the  moment 
we'll  leave  the  choice  of  arbitrary.  The  error  at  the 

end  of  the  r^^  iteration  is  measured  by  the  following: 


{e}  = {AF}  - {AT} 


(30) 


Assume  that  Eq.  (27),  for  t = t]^  , possesses  a unique  solu- 
tion. Then  for  any  ^(^l)J  which  the  iterative  solution 

converges,  it  converges  to  the  correct  solution.  Apparently, 
there  is  some  flexibility  in  the  choice  of  [^(ti)^  , for  a given 
ideal  model,  which  gives  convergent  solutions. 
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The  global  stiffness  matrix  is  selected  conveniently  in  terms 
of  the  element  stiffness  matrices.  If  the  time  increment  is  taken 
short,  the  ideal  model  of  constitutive  behavior  [Eq.  (14)]  is 
approximated  as  follows : 


A(ti) 


m 

(AF) 


m 

+ {t} 


(31) 


where  [A(t]^)J  is  a measure  of  the  stiffness  of  the  material  in 
the  m^^  element  for  t ■ ti  and  {x}  represents  the  change  in 

ZD 

stress  for  fixed  {F}  , e.g.,  thermal  stress,  relaxation,  etc. 

It  follows  Eqs.  (26)^*®  and  (31)  that 


where 


{At} 


in 

[k] 


in 

{Aq} 


/„ 


jn  'T  jn 


d V, 


m 


'm 


[k] 


IQ  rn  m m 

[B]^[A]  [B]  dV^ 


(32) 


(33) 


is  the  element  tangent  stiffness  matrix.  Eqs.  (15),  (26)®,  and 
(32)  lead  to  the  following  representation  of  the  global  stiffness 
matrix: 


M 


(34) 


The  approach  used  in  the  finite  element  solution  is  to  select 
[A(t)]  and  {T(t)}  based  on  the  particular  ideal  model  (i.e., 
constitutive  relations)  being  studied.  Then  the  global  stiffness 
matrix  is  defined  by  Eqs.  (33)  and  (34)  and  the  numerical  solution 
proceeds  as  described  previously.  Prior  to  presenting  the  sample 
calculations  that  were  performed,  the  constitutive  equation  repre- 
sentation for  three  ideal  models  is  discussed. 
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Hookian  Elastic  Model 


The  "Hookian  Elastic"  material  is  represented  by  the  fol- 
lowing constitutive  equation: 


„iY  . pi  / ('fJ  F - z ) 


(35) 


raY66 


where  is  a constant  elastic  stiffness  tensor  and  the 

reference  configuration  is  the  natural  (i.e,,  unstressed)  configu- 
ration. 

Eq.  (35)  is  the  form  appropriate  for  the  present  purposes, 
but  it  may  be  transformed  easily  to  the  following  more  recogniz- 
able form: 


;eY  „ 


66 


(36) 


In  this  report  the  constitutive- behavior  represented  by 
either  Eq.  (35),  or  (36)  is  what  is  meant  by  "Hookian  Elastic. 

Define  as  follows : 

^lY  c = 

^ k ~ k 

It  follows  Eqs . (35),  (36),  and  (37)  that 


(37) 


,iY  e 


(38) 


is  a tensor  representation  of  [A]  which  is  used  for  the 
finite  element  solution  procedure  for  "Hookian  Elastic"  material. 
In  addition,  it  follows  Eq.  (35)  that 


{?}  - 0 


(39) 
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Isotropic  Elastic  Model 

In  this  report  a specific  form  of  Isotropic  elastic  behavior 
Is  selected  for  this  model: 

4,.  - [k  (J  - 1)  - U if]'’  I* 

where  the  reference  configuration  Is  the  natural  (unstressed) 
configuration,  K Is  a scaler  bulk  modulus,  y Is  a scaler 
shear  modulus , and 

J « Det  F 

(41) 

IB  » tr  Fj  g“® 

A clearer  representation  of  the  behavior  Is  represented  by  the 
Cauchy  stress,  Eq.  (40)  Is  equivalent  to  the  following: 


L 


K (J  - 1)  - y 


Ib"!  „1  . y ni 
«J«J  J 


where 


= F^  F^  e 

j g^j 


It  follows  Eq,  (42)  that 


rpf  1 rpk  i 


3K  (J  - 1) 


y IB  ^1\ 


It's  clear  from  Eq.  (43)  that  this  model  represents  an  Isotropic 
elastic  solid  and  that  K and  y are  bulk  and  shear  modulll, 
respectively. 


Define  as  follows : 


e _ 


It  follows  Eqs . (40),  (44)  that 


- K J I 


, T?  « X 


(2J  - 1)  F^  - (J  - 1)  F 


1 -1  Ei~ 

? V ^ ^ 

ak 


* ^ [44  + I - 2r£ra^)_ 


• m 

is  a tensor  representation  of  [A]  which  is  used  for  the 

finite  element  solution  procedure.  In  addition,  it  follows 

Eq.  (40)  that 

{t}  = 0 (4( 


Compressible  Newtonian  Fluid  Model 

In  this  report  the  following  constitutive  equation  represents 
a compressible  Newtonian  fluid: 


p(J)  J + 


2 

I 


+ (pj) 


where  p(J)  is  a scaler  valued  function  and  the  superposed  dot 
represents  the  material  derivative. 
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Eq.  (47)  may  be  represented  In  terms  of  Cauchy's  stress 
tensor  using  kinematical  Identities  and  some  algebra: 

- - p(J)  g^^  + y - I L®  g^J)  (48) 

where  L is  the  velocity  gradient. 

Further  clarification  follows  the  partitioning  of  Eq.  (48) 
into  the  trace  of  T and  the  deviator  of  T : 

T®  - - 3 p(J) 
m 

(49) 

It  follows  Eq.  (49)  that  the  trace  of  T is  a scaler  function  of 
J , the  value  of  p(J)  is  the  pres^ure,  there  is  no  bulk  viscos- 
ity, and  shear  flow  induces  shear  stress  with  a constant  viscosity 
y . 

It's  convenient  to  rearrange  Eq.  (47)  for  numerical  computa- 
tions. After  substituting 


j - J 7 j -// 


Eq.  (47)  may  be  presented  as  follows: 


^Ra 


-li 

- PU)  J 


.1  J"  6 j 


-iBi 


6Vi 


^ak  - l^k  VFb 


(50) 


r 
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In  the  finite  element  representation  it  is  assiamed  that  the 
nodal  coordinates  vary  linearly  with  time  during  an  increment. 

It  follows  Eq,  (17)^  that  F is  constant  during  an  increment. 
Holding  F fixed  in  Eq.  (50),  the  time  dependence  of  is 

seen  to  be  implicit — i.e.,  depends  on  F(t)  — during  an  incre- 
ment . Since  F is  constant  during  an  increment , it  may  be 
expressed  as  follows : 

k ^^6 

"s  - -ST  » 

Let  be  defined  as  follows: 


‘ ak 


-f  y J (gf  't  , - I ^ 

I ^k  j a ak  3 k a 


j A ni 

Again  is  a tensor  representation  of  [A]  and  it 


follows  Eq.  (50)  that 


xn 

(t)  = 0 


2 . Sample  Problems 
Test  Case  #1 

Sample  for  Hookian  Elastic  Material 

Consider  an  orthotropic  body  in  plane  stress  whose  elastic 
modulii  are  given  the  following  values  relative  to  principal 
material  directions . 


10°  psi  V, 


.1  E««  » 10^  psi 


Gt  ~ * .5 X 10^  psi  V,«  “ .01 


In  plane  stress  the  stiffness  tensor  (K)  for  a Hookian 
Elastic  material  may  be  reduced  to  two  dimensional.  The  nonzero 
components  of  K , which  correspond  to  the  values  of  moduli! 
given  by  Eq.  (54),  are 


K 


1111 


10°  psi 


K 


1122 


K 


2211 


10^  psi 


k2222  . 3^q5 

j,1212  . j,1221  . j^2112  ^ j,2121  ^ 5 , 


Now  consider  the  following  problem.  The  body  is  rectangular 
with  unit  width  and  length  three  times  the  width.  It  is  referred 
to  rectangular  Cartesian  coordinates  and  is  in  a natural  state  at 
t * 0 . Loading  (traction  and  boundary  motion)  is  specified  as  a 
function  of  time  (i.e.,  for  t>0).  At  t = 0 the  body  occupies 
the  following  space: 


0 < < 1 


0 5 X2  ^ 3 


The  problem  is  specified  by  the  following  surface  loading; 

on  X^  = 0 and  X^  = 1 , surface  traction  is  zero 
2 

on  X = 3 , nonzero  surface  traction  is  specified: 
sj(t)  * sin  e(t)(l  + 2 X 10-2t)^/2(l03  - l)  t 
S^(t)  = cos  e(t)(l  + 2 X 10-2t)^/2(l03  - l)t 
2 

on  X = 0 , motion  is  specified: 

x^(t)  » cos  e(t)(l  - 2xl0-^t)^^^  X^ 


x^(t)  * - sin  e 


(t)( 


1 - 2 X 10-^t)  X^ 


(56) 


(57) 


Figure  1 illustrates  the  surface  loadings. 
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Figure  1.  Surface  loadings  for  Test  Case  //I 


The  following  motion  is  an  exact  solution  to  the  stated 
problem. 

x^(t)  « cos  e(t)(l  - 2 • lO'^t)^^^ 

+ sin  0(t)(l  + 2 • 

(58) 

x^(t)  - - sin  e(t)(l  - 2 • 10‘^t)^^^  X^ 

+ cos  e(t)(l  + 2 • 10"^t)^''^  X^ 


The  deformation  gradient  and  stress  are  easily  computed  from  the 
motion  given  by  Eq,  (58)  and  the  stiffness  tensor. 

The  problem  represented  by  Eq.  (57)  was  solved  numerically 
using  the  finite  element  technique  described  earlier.  It  was  run 
for  a duration  of  1 sec.  Several  variations  were  run,  each  with 
a constant  rate  of  rotation  i.e.,  represented  by  6(t)  . The 
rates  of  rotation  studied  ranged  from  zero  to  90®/sec.  Also  of 
interest  is  the  tradeoff  between  time  step  (i.e.,  increment  size) 
and  the  number  of  iterations  required  for  convergence.  The  con- 
vergence criterion  used  is  that  the  iterative  change  in  each  nodal 
coordinate  be  less  than  .17d  of  the  maximum  coordinate  change 
during  the  time  step.  Figure  2 presents  the  finite  element  model 
used  in  the  study,  and  Table  1 summarizes  the  convergence  study. 

The  accuracy  obtained  is  excellent.  A comparison  of  the 
numerical  results  and  analytic  solution  for  6 = 30®/sec  is 
presented  by  Figure  3. 

Test  Case  #2 

Generalized  Shear 

Green  and  Atkins®  present  an  analytic  solution  for  "Gener- 
alized Shear"  of  an  isotropic  incompressible  elastic  material, 
defined  by  the  following  constitutive  equation: 
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6 Elements 
16  Degrees  of  Freedom 


Figure  2.  Test  Case  #1  - finite  element  model 
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TABLE  1 
TEST  CASE  n 

NUMERICAL  SOLUTION  TABLE  OF  CONVERGENCE 


ROTATION  RATE 
(deg/sec) 


TIME  STEP 
(sec) 


NUMBER  OF 
ITERATIONS 


(59) 


“-pr  TVf 

where  the  deformation  is  restricted  so  that  Det  F * 1 , and  p 
is  a scaler  which  is  independent  of  straining.  Now  consider  a 
rectangular  two-dimensional  body  subject  to  plane  deformation. 

At  t “ 0 the  body  is  in  a natural  state  and  occupies  the  following 


space : 


0 < < .25 L 

( 

0 5 X^  < L 

On  the  four  boundaries  the  following  loading  is  specified  as  a 
function  of  dimensionless  time: 

2 

on  X “ 0 the  boundary  is  fixed 


on  X" 


s^(t)  - u I 


s2(t) 


on  X 


on  X^  - .25L 


- - I 


Figure  4 illustrates  the  surface  loadings. 
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Figure  4.  Test  Case  /A  2 - surface  loadings 


An  analytic  solution  to  the  problem  given  by  Eq,  (61)  is 
given  by  the  following: 


x^(t)  ^ ^ 
L L 


(62) 


The  problem  represented  by  Eq.  (61)  was  solved  numerically 
using  the  finite  element  solution  technique  described  earlier. 

The  constitutive  equation  presented  by  Eq.  (59)  is  an  approximate 
representation  of  rubbery  materials  and  is  restricted  to  isochoric 
deformations.  For  the  numerical  study  the  constitutive  equation 
presented  by  Eq.  (AO)  was  used.  In  rubbery  material  the  ratio  of 
bulk  modulus  to  shear  modulus  is  of  the  order  of  10^,  therefore 
in  the  numerical  study 

- - 10^  (6A) 

y 
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The  finite  element  model  used  for  this  study  is  presented  by 
Figure  5.  The  problem  was  run  for  a duration  of  t*l.  The 
tradeoff  between  time  step  size  and  the  number  of  iterations 
required  for  convergence  is  sxammarized  by  Table  2.  A 1%  con- 
vergence tolerance  means  the  solution  for  any  time  step  was  iter- 
ated vintil  the  change  in  each  nodal  coordinate  from  an  iteration 
was  less  than  1%  of  the  maximum  value  of  any  nodal  displacement 
during  that  time  step.  The  numerical  results  compared  excellently 
with  the  analytic  solution.  For  example,  the  comparison  of  three 
parameters  is  given  by  Figure  6.  (|)  is  the  angle  illustrated  by 

the  sketch  on  the  figure,  Q is  the  x^  coordinate  of  the  top 
comer  node,  and  is  computed  at  the  centroid  of  the  upper 

most  right  element  on  Figure  5. 


Test  Case  #3 

Pure  Shearing  Flow  of  a Compressible  Newtonian  Fluid 
Consider  the  plane  flow  represented  by 


- a X 


2 2 
v^  = + ax‘‘ 


for  t > 0 


(65) 


v^  = 0 


where  a is  constant  and  v is  the  particle  velocity. 

Such  a flow  is  possible  in  a Newtonian  fluid  and  the  corre- 
sponding stresses  are  computed  using  Eq.  (48): 


T 

T 

T 

T 


11 

22 

33 

12 


- p(J)  - 2ya 

- p(J)  + 2va 

- P(J) 

.j,23  ^ ,j,31  _ Q 


(66) 
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Figure  5.  Test  Case  - finite  element  model. 
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TABLE  2 
TEST  CASE  #2 


NUMERICAL  SOLUTION  TABLE  OF  CONVERGENCE 


TIME  STEP 

CONVERGENCE 

TOLERANCE 

NUMBER  OF 
STEPS 

AVG  ITERATIONS 
PER  STEP 

.001 

] 

[% 

1000 

3 

.01 

- 100 

5 

.02 

1 

50 

7 

.025 

AO 

8 

.0^1 

- 

> 50 

.05 

> 50 

1.0 

1 

BLOWS  UP 
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Now  consider  a rectangular  body  subject  to  plane  deformation. 
Let  the  body  occupy  the  following  space  at  t = 0 ; 


0 < < 1.5 

0 < < 2 


(67) 


with  the  initial  velocity  field  given  by 


at  t = 0 


(68) 


Furthermore  subject  the  boundaries  of  the  body  to  the  following 
mixed  specification  for  t > 0 . 

on  X^  = 0 v^  = 0 * 0 

on  X^  = 1.5  v^(t)  = - 1.5  ae’^^ 

s|  = 0 

(69) 

on  X^  = 0 v^  = 0 sj  = 0 

on  X^  = 2 v^(t)  = 2ae®^ 


Figure  7 illustrates  the  surface  motion  specified. 

It  can  be  shown  that  the  problem  stated  by  Eqs . (67),  (68), 
(69)  has  a solution  given  by  Eq.  (65)  with  stresses  given  by 
Eq.  (66).  It  also  follows  that  the  motion  of  the  body  depends 
only  on  the  motion  specified  at  the  boundaries  and  is  independent 
of  the  viscosity  v . In  terms  of  material  coordinates  the 
motion  is  given  by  the  following: 
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Motion  specified 


Motion  specified 


Pure  shearing  flow 


Figure  7.  Test  Case  #3  - surface  motion. 


(70) 


It  follows  Eq.  (70)  that  the  deformation  gradient  is  given  by 
the  following: 

= e-®^ 

f22  ^ ^at  For  t ^ 0 (71) 

= 0 

The  problem  defined  by  Eqs . (67),  (68),  (69)  was  solved 
numerically  using  the  finite  element  technique  described  earlier. 
It  was  run  for  a duration  of  1.2  sec  with  a=  .5/sec  . The 
finite  element  model  used  is  presented  by  Figure  8.  The  numer- 
ical results  compared  excellently  with  the  analytic  solution. 

To  illustrate,  the  comparison  is  given  for  some  of  the  parameters 
by  Figure  9.  The  values  of  F^^  , f22  are  taken  from  one  of  the 
elements,  and  , 0^2  current  coordinates  of  one  of 

the  interior  nodes . 


Test  Case  # 4 

Growth  of  an  Elliptic  Cavity 

In  1962,  Berg®  presented  the  results  of  an  interesting  study: 
"The  Motion  of  Cracks  in  Plane  Viscous  Deformation."  This  paper 
is  the  basis  for  the  choice  of  this  sample  problem.  The  problem 
is  especially  suited  for  the  present  study  because  it  involves 
following  a deforming  free-surface  in  a fluid  body. 

The  problem  is  illustrated  by  Figure  10.  An  infinite  sheet 
of  material  is  subject  to  pure  shear  stress  at  infinity.  At 
t*0  there  exists  an  elliptically  shaped  cavity.  The  material 
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Figure  8.  Test  Case  #3  - finite  element  model 


Analytical  Solution 


Figure  9.  Test  Case  #3  - numerical  solution. 
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is  characterized  as  an  incompressible  Newtonian  fluid.  The 
stress  at  infinity  is  held  constant  and  the  motion  of  the  cavity 
is  studied.  Berg  shows  that  an  elliptical  cavity  remains  ellip- 
tical always.  Also,  it  follows  Berg's  analysis  that  the  flow  at 
infinity  (for  pure  shear  stress)  preserves  the  area  of  the  infi- 
nite sheet.  Since  the  constitutive  behavior  is  incompressible, 
is  is  concluded  (here)  that  the  area  of  the  cavity  must  be 
preserved.  If  we  let  a(t)  , b(t)  represent  the  semi-axes  of 
the  ellipse  it  follows  that  the  solution  must  satisfy 

a(t)  b(t)  = constant  (72) 

A form  convenient  for  checking  numerical  results  is  obtained  by 
differentiating  Eq.  (72): 


a(t)  _ a(t)  r,-,. 

t(t)  b(t)  ^ ^ 

It's  difficult  to  gleam  more  quantitative  data  from  the  analytical 
solution. 

The  problem  was  studied  numerically  using  the  compressible 
Newtonian  fluid  model.  The  bulk  becavior  was  assvimed  given  by  the 
the  following: 

p(J)  = 10^  psi  (1  - J)  (74) 


Obviously  an  infinite  body  cannot  be  modeled.  The  finite  element 
model  used  is  shown  by  Figure  11.  Note  the  break  in  the  coordi- 
nates and  that  the  advantages  of  symmetry  are  utilized.  To 
simulate  the  motion  at  infinity  the  following  motion  was  specified 
for  each  node  located  on  the  outer  ring  of  the  model: 


.1(0  -.1x1  e-^/2 

^(t)  « + i e 


(75) 
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The  initial  nodal  velocities  were  assumed  zero;  this  is  obvi- 
ously not  the  correct  initial  conditions  for  the  problem  but  the 
correct  conditions  were  not  available.  Some  of  the  results  are 
summarized  by  Figure  12.  The  lower  figure  illustrates  the  cavity 
shape  for  three  values  of  time.  The  upper  figure  illustrates  how 
well  Eq.  (73)  is  satisfied  by  the  niimerical  results;  and 

are  the  nodal  coordinates  whose  values  correspond  to  the 
semi-axes  of  the  elliptical  cavity.  Note  the  inaccuracy  at  very 
early  time  due  to  the  inaccurate  initial  values . The  correlation 
appears  good  for  later  times. 
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III.  RESIDUAL  STRESS  ANALYSIS 


1 . Technical  Overview 

A principal  reason  which  motivates  the  study  of  residual 
stress  is  to  improve  the  quality  of  modem  composite  materials 
and  their  application  to  component  design.  The  mechanical 
characteristics  (strength,  stiffness,  toughness,  etc.)  are 
measurable  properties  which  are  very  sensitive  to  "quality"; 
therefore  they  are  frequently  used  as  indicators , besides  they 
are  of-themselves  of  primary  importance  to  component  design  for 
many  applications.  Hopefully,  further  study  of  the  development 
of  residual  stress  during  processing  will  lead  to  increased 
understanding.  This  knowledge  can  be  applied  to  the  design  of 
processes  which  yield  better  quality  material. 

The  magnitude  of  residual  stress  that  exists  in  a component 
which  is  made  up  of  composite  material  depends  on  physical 
changes  that  occur  in  the  constituents  during  processing  and 
the  geometric  arrangement  of  the  constituents  (i.e.,  laminate 
layup,  three-dimensional  weave,  etc.).  These  two  contributors 
represent  quite  different  phenomena  and  may  be  understood  best 
when  discussed  individually. 

The  physical  changes  that  occur  in  the  constituents  may  be 
quite  different  when  comparing  different  composite  systems,  but 
each  can  be  reduced  to  a single  common  parameter  which  is  essen- 
tial to  residual  stress  development.  The  essential  parameter  is 
the  dimensional  change  experienced  by  the  constituents.  Examples 
are  the  shrinkage  of  a resin  as  a result  of  polymerization,  the 
shrinkage  of  liquid  metal  on  solidification,  growth  of  graphite 
fibers  (due  to  annealing)  during  graphiterization  cycles,  differ- 
ential thermal  expansion  on  cooling  to  room  temperature,  etc. 

The  other  contributor  is  due  to  the  geometry  of  the  constit- 
uents' arrangement  and  their  different  mechanical  properties.  As 


A1 


far  as  predicting  stresses  is  concerned,  it  is  both  convenient 
(from  a computational  viewpoint)  and  enlightening  (from  an 
understanding  viewpoint)  to  perform  different  levels  of  analysis : 
micromechanics  and  macromechanics . In  micromechanics , the 
heterogeneous  nature  of  a material  is  recognized,  and  stresses 
in  individual  constituents  may  be  identified.  On  the  other  hand, 
macromechanics  describes  a material  as  homogeneous  and  uses 
effective  properties  to  represent  its  elastic  behavior.  In  this 
approach,  the  stress  is  interpreted  as  average  (on  an  area  larger 
than  microdimensions) , and  no  information  about  constituent 
stresses  is  obtained. 

It's  useful  to  partition  the  analyses  of  residual  stress 
into  two  types,  which  are  called  here  Predictive  and  Descriptive. 

By  "Predictive  Analysis"  is  meant  a theoretical/numerical  pre- 
diction of  residual  stresses  which  does  not  utilize  residual 
strain  measurements  in  the  component  under  study.  There  are  many 
published  experimental  studies  of  residual  stress;  these  are 
called  here  Descriptive.  Descriptive  analyses  always  involve  an 
experiment  performed  on  a particular  component.  Generally  only 
strains  can  be  measured  and  theoretical/numerical  predictions  are 
required  to  interpret  the  data  in  terms  of  stress.  Clearly,  both 
types  of  analyses  have  shortcomings , but  both  contribute  toward 
improved  understanding.  The  publications  of  Shaffer  and  Levitsky’’® 
represent  an  ambitious  approach  toward  a predictive  analysis. 
Examples  of  predictive  analysis  in  use  for  modern  composite  mate- 
rials are  published  by  Chamis®.  Descriptive  analyses  of  interest 
are  found  in  the  publications  of  Marloff  and  Daniel'®,  Koufopoulos ' ' 
Daniel  and  Liber' and  Lee,  Sacher,  and  Lewis'®. 

It's  impossible  to  perform  a predictive  analysis  with  assuming 
a model  to  represent  the  physical-chemical  processes  at  work  during 
the  fabrication  of  a component  — even  if  the  assumption  is  only 
implied  by  the  retention  and  neglect  of  various  phenomenon.  In 
order  to  be  definitive  the  following  describes  a model  assumed  here 
for  resin  reinforced  composite  materials: 
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1)  The  resin  is  in  a green  stage  when  a component  is  laid 
up,  and  prior  to  curing  the  component  is  stress -free. 

2)  Curing  is  performed  under  controlled  pressure  at 
elevated  temperatures  (frequently  in  the  range  of 
300  to  400  ®F). 

3)  During  the  cure  time  the  resin  polymerizes.  When  the 
component  is  first  heated  to  the  cure  temperature  the 
degree-of-polymerization  (D-O-P)  is  low  and  the  resin 
flows  easily.  With  time  the  D-O-P  advances  and  the 
resin  viscosity  increases,  making  it  gradually  more 
difficult  for  the  resin  to  flow.  Simultaneously  the 
stress-free  volume  of  the  resin  decreases  (or  increases) 
as  the  D-O-P  advances,  and  it  is  assumed  the  dimensions 
of  the  reinforcement  are  stable. at  the  cure  temperature. 
Consequently,  the  resin  shrinkage  induces  stresses  in 
the  component.  Finally,  with  sufficient  time,  the 
polynerization  is  complete  and  the  final  stress  distri- 
bution in  the  component  depends  on  the  history  of 
shrinkage  and  viscosity  in  addition  to  the  geometry  and 
modulii  of  the  reinforcement. 

4)  After  the  completion  of  polymerization  the  component  is 
cooled  to  room  temperature  and  the  pressure  applied  to 
the  component  is  relieved.  During  the  cool-down  the 
tvo  constituents  undergo  unequal  thermal  contraction. 

'It  . resulting  incompatible  dimensional  changes  during 
the  cool-down  induce  further  stress 

5)  The  residual  stress  is  the  stress  remaining  in  the  com- 
ponent at  room  temperature.  It  is  a result  of  incom- 
patible dimensional  changes  which  are  a consequence  of 
both  "polymerization  changes"  and  "cool-down." 

No  one  doubts  the  importance  of  cool-down  for  the  consider- 
ation of  residual  stresses.  It  is  the  only  phenomenon  considered 


in  predictive  analyses  of  modem  composite  materials.  On  the 
contrary  the  effects  of  polymerization  are  ignored  generally  — 
not  through  ignorance  necessarily.  Many  of  the  researchers 
believe  that  creep  or  relaxation  of  the  resin  at  cure  temperature 
is  sufficient  so  that  the  stresses  developed  during  polymer- 
ization are  negligible;  at  least  this  appears  to  be  a common 
opinion  for  epoxy  resin.  Yet  there  appears  to  bte  little  evidence 
in  the  available  descriptive  analyses  to  support  such  an  assxjmp- 
tion. 


The  conclusion  (made  here)  that  the  effects  of  the  polymer- 
ization change  have  not  been  identified  in  descriptive  analyses 
is  not  immediately  obvious  upon  reviewing  the  available  docu- 
mentation. On  first  reflection  one  easily  may  think  that  the 
results  of  Daniel  and  Liber^^  measure  directly  the  effects  of 
the  polymerization  change.  They  implanted  strain  gages  between 
the  lamina  of  a laminate  in  the  green  stage  and  monitored  strain 
output  during  the  polymerization  of  the  resin.  The  following 
idealized  description  of  the  phenomena  sheds  doubt  on  the  useful- 
ness of  the  imbedded  strain  gage  as  a measure  of  polymerization 
changes . 

The  idealized  model  is  illustrated  schematically  by  Figure  13. 
The  sketch  on  the  upper  right  illustrates  that  the  viscosity  of 
the  resin  increases  with  time  during  curing  (i.e.,  implicitly 
through  D-O-P)  and  the  stress-free  volume  decreases  with  time. 

The  sketch  on  the  upper  left  suggests  that  for  a period  (say 
t < t]^)  the  resin  flows  readily;  therefore,  stresses  due  to  the 
incompatibility  are  relaxed  quickly.  Suppose  at  time  t = t^  the 
viscosity  has  grown  sufficiently  so  that  further  development  of 
stress  is  retained,  i.e.,  further  straining  of  the  resin  is  con- 
strained by  the  fibers.  Suppose  the  stress-free  volume  is  V-^ 
at  t * t]^  and  Vp  is  the  value  for  fully  cured  resin.  It 
follows  that  the  final  residual  stress  will  be  proportional  to 
Vp  - Vi  . The  uniaxial  strain  corresponding  to  Vp  - V]^  is  not 
observable  to  the  imbedded  strain  gage;  Vp  - is  a mechanical 
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strain  of  the  resin  which  is  prevented  from  occurring  by  the 
much  stiffer  fibers.  Hence,  the  residual  stress  developed  during 
polymerization  is  due  to  a lack  of  (potential)  straining  which 
would  occur  were  it  not  for  the  constraint  of  the  stiffer  fibers. 
Imbedded  strain  gages  sense  only  the  (total)  strain  that  dees 
occur;  therefore,  such  strain  data  is  insensitive  to  the  residual 
stress  that  develops  during  polymerization 

In  the  absence  of  direct  descriptive  analyses,  other  indi- 
cators of  the  importance,  or  unimportance,  of  polymerization 
changes  may  be  sought.  For  example,  the  total  incompatibility 
may  be  assessed  in  the  absence  of  relaxation.  If  the  total  incom- 
patibility is  unimportant  for  a particular  resin,  it  appears  safe 
to  ignore  it.  Such  an  assessment  is  summarized  by  Table  3. 

Three  resins  were  selected.  Epoxy  and  phenolic  are  two 
resins  for  which  there  is  considerable  processing  experience  and 
interest.  Polyimide  is  selected  because  presently  there  is  con- 
siderable  interest  in  its  application,  e is  the  uniaxial 
thermal  expansion  of  the  resin  due  to  cooling  from  the  deposition 
temperature  to  room  temperature,  is  the  uniaxial  dimensional 

change  due  to  polymerization,  and  is  the  uniaxial  thermal 

expansion  of  the  reinforcement  due  to  cooling  from  the  curing 
temperature  to  room  temperature.  Of  the  properties  used  in 
Table  3,  is  the  least  available.  The  value  used  for  epoxy 

is  obtained  from  Emerson  & Cvuning  Inc.  (Catalog)  for  Stycast  1269A. 
The  value  used  for  phenolic  is  an  estimate  based  on  private  com- 
munication with  Dr.  Paul  Juneau  (General  Electric  Company, 
Switchgear  Equipment  Business  Division) . No  estimate  was  ob- 
tained for  polyimide  even  though  several  researchers  were  con- 
tacted (e.g.,  Mr.  Benson  Dexter,  NASA-LRC;  Dr.  Terry  Saint  Clare, 
NASA-LRC ; Dr.  Tito  Serafini,  NASA-Lewis ; and  Dr.  Joe  Augl , 

NSWC-WOL) . 

The  incompatibility  due  to  cool-down  is  proportional  to  the 
differential  thermal  expansion  and  that  due  to 


1 
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TABLE  3.  ASSESSMENT  OF  POLYMERIZATION  DIMENSIONAL  CHANGE 


47 


POLYIMIDE  - UNKNOWN 


polymerization  is  represented  by  . The  ratio  of  to  AeT 

(shovm  in  the  last  column)  is  a measure  of  the  importance  of 
polymerization  relative  to  the  cool-down.  For  the  reinforcements 
considered  the  polymerization  changes  for  epoxy  resin  are  in 
excess  of  60%  of  that  due  to  cool-down  and  for  phenolic  resin  are 
in  excess  of  400%  of  that  due  to  cool-down.  It  follows  that  the 
contribution  of  pol3mierization  changes  to  the  residual  stress 
cannot  be  shown  negligible  on  the  basis  of  a conservative  estimate 
such  as  siommarized  by  Table  3. 

2 . Predictive  Analysis 

This  section  summarizes  an  initial  step  toward  the  develop- 
ment of  a predictive  analysis  for  resin  reinforced  composite 
materials . The  approach  taken  is  to  develop  a micromechanical 
analysis  (i.e.,  the  resin  and  reinforcement  are  distinguished  in 
the  analysis).  It’s  assumed  that  the  phenomena  which  occur 
during  polymerization  are  "slow','"  so  that  the  approximations  of 
creeping  flow  are  appropriate;  in  this  case  the  equations-of- 
motion  reduce  to  equilibrium  equations.  The  physical-chemical 
processes  that  occur  during  the  curing  of  the  resin  — and  are 
responsible  for  the  development  of  stress  during  that  time  — are 
reflected  by  the  constitutive  equations  which  are  selected  to 
represent  the  resin  behavior.  The  finite  element  technique, 
which  is  described  in  the  first  section  of  this  report,  is  used 
to  perform  numerical  studies.  The  other  field  equations,  besides 
the  constitutive  equation,  are  contained  inherently  in  the  finite 
element  technique  and  need  no  special  attention  here. 

The  following  list  is  a goal  for  the  phenomena  to  be  repre- 
sented by  the  constitutive  equations: 

1)  The  reinforcement  is  dimensionally  stable  at  the  cure 
temperature  and  behaves  elastically. 

2)  The  resin  starts  liquid-like  with  low  viscosity;  the 
viscosity  increases  with  advance  of  D-O-P  until  at  the 
end  of  cure  the  resin  behaves  solid-like. 
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3)  The  stress-free  volvnne  of  the  resin  varies  continnously 
with  the  D-O-P. 

It  is  useful  to  represent  time -dependent  models  of  mechanical 
behavior  by  springs  and  dash-pots.  The  following  sketch  illus- 
trates pictorially  a representation  of  the  phenomenon  described 
above : 


Bulk  Behavior 


Shear  Behavior 


tF  Is 


K(D-O-P) 


JS  fF 

ill 


M (D-O-P) 


G(D-O-P) 


U (D-O-P) 
F = K(8-U) 


M dt  \ G / 


In  the  representation  of  bulk  behavior  the  bulk  modules  (K) 
may  depend  on  the  D-O-P,  and  the  base  movement  (-U)  represents 
the  change  of  stress -free  volume  with  D-O-P.  In  shear  the 
transition  from  fluid-like  to  solid-like  behavior  is  controlled 
primarily  by  the  viscosity  y . For  small  y the  model  will 
creep  readily  under  applied  stress,  and  the  elastic  recovery  will 
be  very  small.  For  very  large  y an  increment  of  strain  is 
accommodated  by  stretch  of  the  spring,  which  can  be  recovered 
elastically.  Actually,  the  magnitude  of  y , small  or  large, 
must  be  measured  relative  to  the  magnitude  of  G and  a measure 
of  time. 


The  choice  of  a multiaxial  invariant  constitutive  equation, 
for  which  the  spring  and  dash-pot  model  above  is  an  analogue,  is 
not  unique.  Following  Eq,  (43),  B is  selected  as  the  measure 
of  strain  and  the  corotational  derivative  is  selected  for  the 
time  flux;  the  constitutive  equation  is  represented  as  follows: 


rt 

- 3K(t)(j(t)  - Jjj(t)) 

Ti,(0 

u(t) 

Tij(t)‘ 

G(t) 

(76) 

where  t is  the  time,  t*0  corresponds  to  the  start  of  cure, 
the  reference  configuration  corresponds  to  t * 0 (it  follovrs 
that  (0)  * g^j  , J(0)“1  ),  the  time  dependence  of  v and 

G is  defined  implicitly  via  the  time  dependence  of  D-O-P,  and 

T^(t) 

Tj^(t)  - Tj(t)  - -5L—  gj  (77) 

and  the  corotational  derivative  of  a second-order  tensor  (A^j) 
is  given  by  the  following: 

if  *ij  - ^ - (“Ki  ''Kj ) (’8) 

d * 

where  is  the  material  derivative  and  w^  is  the  spin 

tensor . 

It's  convenient  to  define  a tensor  Aj^j  (t)  which  is  a 
measure  of  the  integrated  viscous  strain: 

T ' 

Atj(t)  . (79) 

t 

It  follows  Eq.  (79)  that  A^j  is  symmetric  and  trAj^j(t)«0  . 


The  constitutive  equations  may  now  be  represented  as  follows: 


) - 3K(t)(j(t)  - 

T'  (t)  r ~1 


(80) 


For  the  finite  element  technique  the  material  derivatives 
must  be  expressed  in  terms  of  increments.  The  following  repre- 
sentations are  used: 


dt 

^(tn) 


(^n) 

At 


(81) 


-1. 


Eqs . (79),  (80),  and  (81)  complete  the  formulation  of  the 
constitutive  equation.  The  global  stiffness  matrix  required  for 
the  iterations  is  generated  by  Eq.  (45). 

The  following  idealized  cure  process  is  assumed  in  this 
s tudy : 

1)  The  component  is  cured  for  2 hours  at  350®F  and  a pres- 
sure of  85  psia.  Time  is  measured  from  the  start  of 
the  cure  cycle,  and  the  stresses  throughout  the  compo 
nent  are  equivalent  to  the  applied  hydrostatic  pressure 
at  t “ 0 . 

2)  The  environmental  pressure  is  reduced  to  atmospheric 
pressure  in  the  next  10  minutes  and  then  the  comnonent 
is  cooled  to  75®F  during  the  next  2 hours. 
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The  following  tabulation  summarizes  the  material  functions 
used  for  the  numerical  study,  and  it  is  assumed  that  they  vary 
linearly  with  time  between  the  instants  given: 

EPOXY 

TIME 


0 

2 hrs 

2 hrs  - 10  min 

4 hrs  - 10  min 

K(t)(10^  psi) 

.16 

.16 

.16 

.75 

G(t)(10^  psi) 

.076 

.076 

.076 

.18 

y(t)(10^  psi-min) 

1.1 

1.1  X 10^ 

1.1 X 10^ 

1.7  X 10  ^ 

JN(t)  - 1 

.4375  X 10-3 

-.01757 

-.01757 

-.04017 

PHENOLIC 

TIME 


0 

2 hrs 

2 hrs  - 10  min 

4 hrs  - 10  min 

K(t)(10^  psi) 

.16 

.16 

.16 

.75 

G(t)(10^  psi) 

.076 

.076 

.076 

.18 

y(t)(10^  psi-min) 

1.1 

1.1 X 10^ 

1.1 X 10^ 

1.7  X 10^ 

J^Ct)  - 1 

.4375  X 10-3 

-.17964 

-.17964 

-.2141 

Some  rationale  for  the  values  chosen  is  in  order.  Consider 
the  epoxy  characterization.  Wang^®  published  mechanical  test  data 
for  cured  epoxy.  The  values  of  K and  G at  4 hrs  - 10  min  were 
selected  to  correspond  to  Wang's  tensile  test  data  at  72®F;  the 
values  of  K and  G at  2 hrs  - 10  min  and  at  2 hrs  were  selected 
to  correspond  to  his  tensile  test  data  at  375®F.  The  value  of  y 
at  4 hrs  - 10  min  was  selected  to  correspond  to  Wang's  tensile  creep 
data  at  1000  psi  and  375®F;  the  value  of  y at  2 hrs  - 10  min  and 
2 hrs  was  selected  to  correspond  to  his  tensile  creep  data  at 
2000  psi  and  73®F.  The  values  of  K , G , and  y at  0 hrs  were 
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assumed  with  no  experimental  justification;  the  value  of  y for 
the  fluid-like  vincured  state  of  the  material  was  chosen  very 
small  relative  to  its  value  for  the  cured  state.  Intuitively 
it  appears  that  what  happens  during  the  cure  cycle  is  less  sensi- 
tive to  K and  G than  to  y , and,  since  no  data  were  avail- 
able for  these  properties,  the  values  of  K and  G were  selected 
equal  to  their  fully  cured  values . The  value  of  Jjj  at  0 hrs 
was  selected  so  that  the  stress  in  the  resin  corresponds  to  the 
applied  hydrostatic  pressure;  the  value  of  Jjj  at  2 hrs  and 
2 hrs  - 10  min  was  selected  to  correspond  to  the  polymerization 
dimensional  change  given  by  Table  3.  The  value  of  Jjj  at  4 hrs  - 
10  min  was  selected  to  correspond  to  the  differential  thermal 
expansion  (Ae*^)  assuming  graphite  reinforcement  given  by  Table  3. 

Mechanical  test  data  similar  to  Wang's  for  epoxy  were  not 
available  for  phenolic.  The  values  of  K, , G , and  y were 
selected  to  match  those  used  for  epoxy.  The  values  of  Jjj  were 
selected  using  the  same  rationale  which  was  used  for  the  epoxy. 

The  finite  element  model  of  the  material  is  represented  by 
Figure  14.  Attention  is  concentrated  on  the  resin  surrounding  a 
single  fiber  in  a laminate  of  50  volume  %.  It's  assigned  that  the 
model  represents  a fiber  in  an  interior  lamina  of  a multi - 
directionally  reinforced  laminate.  It's  also  assumed  that  the 
fibers  are  rigid  relative  to  the  resin  (e.g..  Young's  modulus 
of  a graphite  fiber  is  100  times  that  of  epoxy) . Advantages  of 
symmetry  are  utilized  by  modeling  only  1/4  of  the  fiber  neighbor- 
hood. A thin  plate  of  stiff  (Hookian  Elastic)  material  (repre- 
sented by  elements  1 through  6)  is  added  to  the  top  edge  of  the 
resin;  it's  used  to  apply  and  remove  the  controlled  hydrostatic 
pressure  and  maintain  flatness  of  the  top  edge.  The  circular 
boundary  (contiguous  with  the  fiber)  is  held  fixed.  On  the  other 
boundary  the  normal  displacement  and  the  shear  traction  are  set 
to  zero. 
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The  results  are  summarized  by  Figures  15  and  16.  A review 
of  the  printout  shows  that  the  stresses  are  maximum  generally  in 
the  neighborhood  of  Element  #40  and  that  the  state  of  stress  is 
well  represented  by  the  mean  stress  (l.e.,  T™/3  ) — because  the 

deviatoric  component  is  but  a few  percent  of  the  mean  stress . 

The  cases  were  performed  using  an  increment  size  of  one  minute 
and  a convergence  criterion  of  17,. 

The  time  history  of  the  mean  stress  in  Element  #40  is  pre- 
sented for  the  epoxy  resin  by  Figure  15  and  for  the  phenolic  resin 
by  Figure  16.  In  order  to  present  clearly  the  effects  of  polymer- 
ization dimensional  change,  two  cases  were  studied  for  each  resin; 
one  case  includes  both  polymerization  and  cool-down,  and  the 
second  case  includes  only  the  cool-down.  The  results  Indicate  a 
considerable  effect  on  the  residual  stress  — even  for  the  epoxy 
resin. 

The  magnitude  of  the  predicted  stress  is  very  large  — approx- 
imately 507o  greater  than  the  reported  laniaxial  tensile  strength 
of  epoxy  for  the  case  without  polymerization  and  the  predicted 
stress  is  relatively  greater  for  the  other  cases. 

The  numerical  results  presented  by  Figures  15  and  16  are  the 
initial  calculations  performed.  Obviously,  the  predicted  values 
are  sensitive  particularly  to  the  constitutive  equations  and  the 
values  of  the  material  parameters  used  to  represent  the  material 
behavior.  More  numerical  studies  and  correlation  with  experi- 
mental data  are  required  to  establish  the  credibility  and  utility 
of  the  technology.  Even  so,  one  is  tempted  to  derive  conclusions 
from  the  numerical  results  even  at  this  early  time  in  the  develop- 
ment of  the  technology.  In  that  case  it's  unavoidable  but  to  con- 
conclude — for  the  composites  studied — the  following: 

1)  One  may  expect  damage  in  the  laminate  — in  the  form  of 
voids  or  cracks  in  the  resin  and/or  debonding  from  the 
fiber — located  between  parallel  fibers  in  a lamina. 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 
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Figure  15.  Mean  stress  in  Element  y/40  for  epoxy  resin 
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Figure  16.  Mean  stress  in  Element  #40  for  phenolic  resin. 
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2)  The  dimensional  changes  associated  with  polymerization 
of  the  resin  should  not  be  neglected  in  the  analysis  of 
residual  stress  or  damage  resulting  therefrom. 
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A comparison  of  numerical  predictions  with  some  exact  analytic  solutions  is  presented, 
and  the  comparison  is  excellent.  The  second  study  concerns  residual  stress.  Its 
development  in  polymeric  filamentary  reinforced  composite  materials  is  discussed  in  terms 
of  differential  dimensional  changes  that  occur  during  processing.  Two  phenomena  which 
contribute  to  the  development  of  residual  stress  are  considered:  dimensional  change  of 
the  resin  due  to  polymerization  and  differential  thermal  expansion  of  the  resin  and 
reinforcement  upon  cooling  from  the  cure  temperature  to  room  temperature.  Results  of 
the  study  to  date  do  not  indicate  that  dimensional  changes  due  to  polymerization  may  be 
neglected  in  residual  stress  assessments.^ 
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